Method and apparatus for processing medical images

ABSTRACT

A method of processing medical images, performed by a computer processor and comprising the steps of: (a) obtaining one or more atlases containing one or more images in which one or more anatomical features have been labelled with label data; (b) obtaining a plurality of unlabelled images; (c) comparing the labelled and unlabelled images and selecting one or more unlabelled images that most closely resemble(s) one or more of the labelled images; (d) to each of those selected image(s), propagating label data from one or more of the closest of the labelled images, thereby labelling the corresponding anatomical feature(s) of each of the selected image(s) and causing the selected image(s) to become labelled image(s); and (e) iteratively repeating from step (c), thereby labelling others of the unlabelled images.

This invention relates to a method and corresponding apparatus forprocessing medical images. It is particularly suitable, but by no meanslimited, for processing magnetic resonance images, for example of thehuman brain.

BACKGROUND TO THE INVENTION

The automated extraction of features from magnetic resonance images(MRI) of the brain is an increasingly important process in neuroimaging.Since the brain anatomy varies significantly across subjects and canundergo significant change, either during aging or through diseaseprogression, finding an appropriate way of dealing with anatomicaldifferences during feature extraction has gained increasing attention inrecent years.

Amongst the most popular methods for dealing with this variability areatlas-based approaches. In the context of the present work, an “atlas”is a dataset (which may be a 3D image, a 2D image, images of anydimension, or a set of images) having annotations or labels in order toidentify points, regions or structures within the image.

Atlas-based approaches assume that the atlases can encode the anatomicalvariability either in a probabilistic or statistical fashion. Whenbuilding representative atlases, it is important to register all imagesto a template that is unbiased towards any particular subgroup of thepopulation. Two approaches using the large deformation diffeomorphicsetting for shape averaging and atlas construction have been proposed byAvants and Gee (2004) and Joshi et al. (2004). Template-free methods forco-registering images form an established framework for spatial imagenormalization. In a departure from approaches that seek a singlerepresentative average atlas, two more recent methods describe ways ofidentifying the modes of different populations in an image dataset(Blezek and Miller, 2007; Sabuncu et al., 2008).

To design variable atlases dependent on subject information, a varietyof approaches have been applied in recent years to the problem ofcharacterizing anatomical changes in brain shape over time and duringdisease progression. Davis et al. (2007) describe a method forpopulation shape regression in which kernel regression is adapted to themanifold of diffeomorphisms and is used to obtain an age-dependentatlas. Ericsson et al. (2008) propose a method for the construction of apatient-specific atlas where different average brain atlases are builtin a small deformation setting according to meta-information such assex, age, or clinical factors.

Methods for extracting features or biomarkers from magnetic resonance(MR) brain image data often begin by automatically segmenting regions ofinterest. A very popular segmentation method is to use label propagationwhich transforms labels from an atlas image to an unseen target image bybringing both images into alignment. Atlases are typically, but notnecessarily, manually labelled. Early work using this approach wasproposed by Bajcsy et al. (1983) as well as more recently Gee et al.(1993) and Collins et al. (1995). The accuracy of label propagationstrongly depends on the accuracy of the underlying image alignment. Toovercome the reliance on a single segmentation, Warfield et al. (2004)proposed STAPLE, a method that computes for a collection ofsegmentations a probabilistic estimate of the true segmentation.Rohlfing et al. (2004) demonstrated the improved robustness and accuracyof a multi-classifier framework where the labels propagated frommultiple atlases are combined in a classifier fusion step to obtain afinal segmentation of the target image. Label propagation in combinationwith classifier fusion was successfully used to segment a large numberof structures in brain MR images by Heckemann et al. (2006).

Due to the wide range of anatomical variation, the selection of atlasesbecomes an important issue in multi-atlas segmentation. The selection ofsuitable atlases for a given target helps to ensure that theatlas-target registrations and the subsequent segmentation are asaccurate as possible. Wu et al. (2007) describe different methods forimproving segmentation results in the single atlas case by incorporatingatlas selection. Aljabar et al. (2009) investigate different similaritymeasures for optimal atlas selection during multi-atlas segmentation.Van Rikxoort et al. (2008) propose a method where atlas combination iscarried out separately in different sub-windows of an image until aconvergence criterion is met. These approaches show that it ismeaningful to select suitable atlases for each target imageindividually. Although an increasing number of MR brain images areavailable, the generation of high-quality manual atlases is alabour-intensive and expensive task (see e.g. Hammers et al. (2003)).This means that atlases are often relatively limited in number and, inmost cases, restricted to a particular population (e.g. young, healthysubjects). This can limit the applicability of the atlas database evenif a selection approach is used. To overcome this, Tang et al. (2009)seek to produce a variety of atlas images by utilizing a PCA model ofdeformations learned from transformations between a single templateimage and training images. Potential atlases are generated bytransforming the initial template with a number of transformationssampled from the model. The assumption is that, by finding a suitableatlas for an unseen image, a fast and accurate registration to thistemplate may be readily obtained. Test data with a greater level ofvariation than the training data would, however, represent a significantchallenge to this approach. Additionally, the use of a highly variabletraining dataset may lead to an unrepresentative PCA model as thelikelihood of registration errors between the diverse images and thesingle template is increased. This restriction makes this approach onlyapplicable in cases where a good registration from all training imagesto the single initial template can be easily obtained.

Atlas-based segmentation benefits from the selection of atlases similarto the target image (Wu et al., 2007; Aljabar et al., 2009). However, inpractice, the initial atlases may only represent a specific subgroup ofthe target image population.

There is therefore a desire to be able to propagate a relatively smallnumber of atlases through to a large and diverse set of MR brain imagesexhibiting a significant amount of anatomical variability.

Prior work where automatically labelled brain images were used to labelunseen images did not result in an improvement of segmentation accuracyover direct multi-atlas propagation. In (Heckemann et al., 2006), whenmultiple relatively homogenous atlases were propagated to randomlyselected intermediate images that were used as single atlases for thesegmentation of unseen images, the resulting average Dice overlaps withmanual delineations were 0:80, compared with 0:84 for direct multi-atlaspropagation and fusion. In a second experiment, single atlases werepropagated to randomly selected intermediate subjects that were thenfurther used for multi-atlas segmentation, resulting in Dice overlapswith manual delineations of 0:78 at best.

Further background art is provided by US 2007/0053589 A1, US2008/0154118 A1 and WO 2009/093146 A1, all of which disclose methods forsegmenting image data.

SUMMARY OF THE INVENTION

According to a first aspect of the present invention there is provided amethod as defined in Claim 1 of the appended claims. Thus there isprovided a method of processing medical images, performed by a computerprocessor and comprising the steps of: (a) obtaining one or more atlasescontaining one or more images in which one or more anatomical featureshave been labelled with label data; (b) obtaining a plurality ofunlabelled images; (c) comparing the labelled and unlabelled images andselecting one or more unlabelled images that most closely resemble(s)one or more of the labelled images; (d) to each of those selectedimage(s), propagating label data from one or more of the closest of thelabelled images, thereby labelling the corresponding anatomicalfeature(s) of each of the selected image(s) and causing the selectedimage(s) to become labelled image(s); and (e) iteratively repeating fromstep (c), thereby labelling others of the unlabelled images.

The term “labelled” should be interpreted broadly, to encompass any kindof delineation, segmentation or annotation of an anatomical feature.Similarly, the term “label data” should be interpreted broadly, toencompass any kind of coding that enables an anatomical feature to bedelineated, segmented or annotated on a medical image.

By virtue of the iterative propagation of label data from the closestlabelled images to the unlabelled images, each unlabelled image can besegmented using structurally-similar atlases. As a consequence,relatively large differences between a labelled image and an unlabelledimage may be broken down into a number of small differences betweencomparatively similar initially-unlabelled images through which thelabel data is propagated, enabling registration errors to be reduced.

Preferable, optional, features are defined in the dependent claims.

Thus, preferably the step of comparing the labelled and unlabelledimages comprises embedding the images into a low-dimensional coordinatesystem. This enables the labelled and unlabelled images to be comparedand the differences to be quantitatively evaluated in acomputationally-efficient manner. In certain embodiments thelow-dimensional coordinate system may be a two-dimensional coordinatespace, thus further simplifying the analysis and processing of thedifferences between the images.

Preferably the step of comparing the labelled and unlabelled imagescomprises defining a set of pairwise measures of similarity by comparingone or more respective anatomical features for each pair of images inthe set of images. Particularly preferably this step further comprisesperforming a spectral analysis operation on the pairwise measures ofsimilarity, although those skilled in the art will appreciate that thereare other ways in which this may be accomplished.

The pairwise measures of similarity may represent the intensitysimilarity between a pair of images, and/or the amount of deformationbetween a pair of images.

Preferably the step of propagating label data comprises propagatinglabel data from a plurality of the closest of the labelled images, basedon a classifier fusion technique. This enables the selected image(s) tobe labelled with greater accuracy.

Preferably the method further comprises, after step (d) and before step(e), a step of performing an intensity-based refinement operation on thenewly-propagated label data, in order to further minimize theaccumulation of registration errors during the labelling process.

The images may be of different subjects. Alternatively, at least some ofthe images may be of the same subject but taken at different points intime, thereby enabling intra-subject variance to be identified andstudied.

The images may be magnetic resonance images, or other medical imagesfamiliar to those skilled in the art.

The method may further comprise labelling an anatomical featurerepresentative of the presence or absence of a condition and using thatfeature to derive a biomarker for that condition. On the basis of thebiomarker, the method may further comprise allocating a subject to adiagnostic category, and/or quantifying a subject's response totreatment, and/or selecting a subject's treatment.

According to a second aspect of the present invention there is providedimaging apparatus arranged to implement a method in accordance with thefirst embodiment of the invention. The imaging apparatus may be amedical scanner, such as an MRI scanner, or some other type.

According to a third aspect of the present invention there is providedimage processing apparatus arranged to implement a method in accordancewith the first embodiment of the invention.

According to a fourth aspect of the present invention there is provideda computer system arranged to implement a method in accordance with thefirst embodiment of the invention.

According to a fifth aspect of the present invention there is provided acomputer program comprising coded instructions for implementing a methodin accordance with the first embodiment of the invention.

According to a sixth aspect of the present invention there is providedcomputer-readable medium or physical carrier signal encoding a computerprogram in accordance with the fifth embodiment of the invention.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the invention will now be described, by way of exampleonly, and with reference to the drawings in which:

FIG. 1 illustrates the process of atlas propagation using our newmethod;

FIG. 2 illustrates results showing the discrimination ability fordifferent chosen feature dimensions among four subject groups (healthyyoung, elderly controls, MCI, AD);

FIG. 3 illustrates the MNI152 brain atlas showing the region of interestaround the hippocampus that was used for the evaluation of pairwiseimage similarities;

FIG. 4 illustrates coordinate embedding of 30 atlases based on healthysubjects and 796 images from elderly dementia patients and age-matchedcontrol subjects;

FIG. 5 illustrates a comparison of segmentation results for the righthippocampus on a transverse slice;

FIG. 6 illustrates the development of segmentation accuracy withincreasing distance from the original set of atlases, with each subsetof images used for evaluation being represented by one bar plot;

FIG. 7 illustrates average hippocampal volumes for manual and automaticsegmentation; and

FIG. 8 is a Bland-Altman plot showing the agreement between volumemeasurement based on manual and automatic segmentation of thehippocampus.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

The present embodiments represent the best ways known to the applicantsof putting the invention into practice. However, they are not the onlyways in which this can be achieved.

Primarily, the present embodiments take the form of a method oralgorithm for processing medical (or other) images. The method oralgorithm may be incorporated in a computer program or a set ofinstruction code capable of being executed by a computer processor. Thecomputer processor may be that of a conventional (sufficiently highperformance) computer, or some other image processing apparatus orcomputer system. Alternatively, the computer processor may beincorporated in, or in communication with, a piece of medical imagingequipment such as an MRI scanner.

The computer program or set of instruction code may be supplied on acomputer-readable medium or data carrier such as a CD-ROM, DVD or solidstate memory device. Alternatively, it may be downloadable as a digitalsignal from a connected computer, or over a local area network or a widearea network such as the Internet. As a further alternative, thecomputer program or set of instruction code may be hard-coded in thecomputer processor (or memory associated therewith) arranged to executeit.

Initial Overview

Our method begins with obtaining one or more pre-existing atlases, inwhich a set of digital images have already been labelled or annotated. Aset of images onto which the labels or annotations are to be propagatedare also obtained, for example from an MRI scanner or another piece ofmedical imaging equipment. The images in question may be of the brain.Alternatively they may be of other parts of the human (or animal) body,such as the knee—for example in order to diagnose osteoarthritis.

The atlas propagation and segmentation process using our new method isdepicted in FIG. 1, which shows five steps. Firstly, in step (1), allthe labelled images (i.e. atlases) and unlabelled images are embeddedinto a low-dimensional manifold. In step (2), the N closest unlabelledimages to the labelled images are selected for segmentation. Then, instep (3), the M closest labelled images are registered to each of theselected images (an example for one selected image is illustrated). Instep (4), intensity refinement is used to obtain label maps for each ofthe selected images. Then, in step (5), steps (2)-(4) are iterated untilfurther images (and preferably all of them) are labelled.

As mentioned earlier, atlas-based segmentation benefits from theselection of atlases similar to the target image. Our method provides aframework where this is ensured by first embedding all images in a lowdimensional coordinate system that provides a distance metric betweenimages and allows neighbourhoods of images to be identified. In themanifold learned from coordinate system embedding, a propagationframework can be identified and labelled atlases can be propagated in astep-wise fashion, starting with the initial atlases, until the wholepopulation is segmented. Each image is segmented using atlases that arewithin its neighbourhood, meaning that deformations between dissimilarimages are broken down to several small deformations betweencomparatively similar images and registration errors are reduced. Tofurther minimize an accumulation of registration errors, anintensity-based refinement of the segmentation is done after each labelpropagation step. Once segmented, an image can in turn be used as anatlas in subsequent segmentation steps. After all images in thepopulation are segmented, they represent a large atlas database fromwhich suitable subsets can be selected for the segmentation of unseenimages. The coordinate system into which the images are embedded isobtained by applying a spectral analysis step to their pairwisesimilarities. As labelled atlases are propagated and fused for aparticular target image, the information they provide is combined with amodel based on the target image intensities to generate the finalsegmentation.

Thus, to propagate an initial set of atlases through a dataset of imageswith a high level of inter-subject variance, a manifold representationof the dataset is learned where images within a local neighbourhood aresimilar to each other. The manifold is represented by a coordinateembedding of all images. This embedding is obtained by applying aspectral analysis step to the complete graph in which each vertexrepresents an image and all pairwise similarities between images areused to define the edge weights in the graph. Pairwise similarities canbe measured as the intensity similarity between the images or the amountof deformation between the images or as a combination of the two.

In successive steps, atlases are propagated within the newly definedcoordinate system. In the first step, the initial set of atlases arepropagated to a number of images in their local neighbourhood and usedto label them. Images labelled in this way become atlases themselves andare, in subsequent steps, further propagated throughout the wholedataset. In this way, each image is labelled using a number of atlasesin its close vicinity which has the benefit of decreasing registrationerror.

In an extension of this technique, one or more scans obtained from thesame subject but at different times (so-called “longitudinal” scans) maybe labelled.

After propagating multiple atlases to each baseline scan, spatial priorsobtained from the multiple atlases may be used to segment not only thebaseline scans (as done initially) but also the longitudinal scans.Hence, this extended technique enables the simultaneous segmentation ofdifferent time points (e.g. day 0, day 3, day 15, etc.), which in turnallows a measurement of the differences between time points.

Thus, images of a subject taken at subsequent time points from thebaseline images can be segmented simultaneously and used to identifyintra-subject variance (i.e. differences in anatomical structure withina single subject but at different time points).

Graph Construction and Manifold Embedding

In order to determine the intermediate atlas propagation steps, allimages are embedded in a manifold represented by a coordinate systemwhich is obtained by applying a spectral analysis step. Spectralanalytic techniques have the advantage of generating feature coordinatesbased on measures of pairwise similarity between data items such asimages. This is in contrast to methods that require distance metricsbetween data items such as multidimensional scaling (MDS). After aspectral analysis step, the distance between two images in the learnedcoordinate system is dependent not only upon the original pairwisesimilarity between them but also upon all the pairwise similarities eachimage has with the remainder of the population. This makes the distancesin the coordinate system embedding a more robust measure of proximitythan individual pairwise measures of similarity which can be susceptibleto noise. A good introduction to spectral analytic methods can be foundin von Luxburg (2007) and further details are available in Chung (1997).

The spectral analysis step is applied to the complete, weighted andundirected graph G=(V, E) with each image in the dataset beingrepresented by one vertex v_(i). The non-negative weights w_(ij) twovertices v_(i) and v_(j) are defined by the similarity s_(ij) therespective images. In the present work intensity based similarities areused. A weights matrix W for G is obtained by collecting the edgeweights w_(ij)=s_(ij) every image pair and a diagonal matrix T containsthe degree sums for each vertex

$d_{ii} = {\sum\limits_{j}{w_{ij}.}}$

The dimension of the feature data derived from a spectral analysis stepcan be chosen by the user. In our work, we tested each dimension for thefeature data in turn and assessed the ability to discriminate betweenthe four subject groups (young, AD, MCI and older control subjects). Thediscrimination ability was measured using the average inter-clusterdistance based on the centroids of each cluster for each featuredimension. For the groups studied, it was maximal when usingtwo-dimensional features and reduced thereafter (see FIG. 2). Wetherefore chose to use the 2D spectral features as a coordinate space inwhich to embed the data.

Image Similarities

In the preferred embodiment of our method, we use an intensity-basedsimilarity between a pair of images I_(i) and I_(j). This similarity isbased on normalized mutual information (NMI) (Studholme et al., 1999)which is with the entropy H(I) of an image I and the joint entropyH(I_(i); I_(j)) of two images defined as

${N\; M\; I_{ij}} = {\frac{{H( I_{i} )} + {H( I_{j} )}}{H( {I_{i},I_{j}} )}.}$

For example, when segmenting the hippocampus, we compute the similaritymeasure between a pair of images as the NMI over a region of interest(ROI) around the hippocampus. The framework is, however, general and auser can choose the similarity measure and region of interestappropriate to the region or structure being segmented. To define theROI, all training images are automatically segmented using standardmulti-atlas segmentation (Heckemann et al., 2006). The resultinghippocampal labels are then aligned to a known brain atlas (e.g. theMNI152-brain T1 atlas (Mazziotta et al., 1995)) using a coarse non-rigidregistration modelled by free-form deformations (FFDs) with a 10 mmB-spline control point spacing (Rueckert et al., 1999) between thecorresponding image and the atlas. The hippocampal ROI are then definedthrough the dilation of the region defined by all voxels which arelabelled as hippocampus by at least 2% of the segmentations. To evaluatethe pairwise similarities, all images are aligned to the known atlasusing the same registrations used for the mask building. FIG. 3 showsthe ROI around the hippocampus superimposed on the brain atlas used forimage normalization.

Segmentation Propagation in the Learned Manifold

In order to propagate the atlas segmentations through the dataset usingthe learned manifold, all images I ε I are separated into two groups,containing the labelled and unlabeled images. These groups are indexedby the sets L and U respectively. Initially, L represents the initialatlas images and U represents all other images. Let d(I_(i); I_(j))represent the Euclidean distance between images I_(i) and I_(j) in themanifold. The average distance from an unlabeled image I_(u) to alllabelled images is:

${\overset{\_}{d}( {I_{u},L} )} = {\frac{1}{L}{\sum\limits_{l \in L}{{d( {I_{u},I_{l}} )}.}}}$

At each iteration, the images I_(u), u ε U with the N smallest averagedistances d(I_(u)) are chosen as targets for propagation. For each ofthese images, the M closest images drawn from I_(l), l ε L are selectedas atlases to be propagated. Subsequently, the index sets U and L areupdated to indicate that the target images in the current iteration havebeen labelled. Stepwise propagation is performed in this way until allimages in the dataset are labelled.

N is an important parameter as it determines the number of imageslabelled during each iteration and therefore it strongly affects theexpected number of intermediate steps that are taken before a targetimage is segmented. M defines the number of atlas images used for eachapplication of multiatlas segmentation. A natural choice is to set M tothe number of initial atlases. Independent of the choice of N, thenumber of registrations needed to segment K images is M×K. The processof segmentation propagation in the learned manifold is summarized inAlgorithm 1:

Algorithm 1: Segmentation propagation in the learned manifold Set L torepresent the initial set of atlases Set U to represent all remainingimages while |U| > 0 do for all I_(u) ∈ U do calculate d(I_(u),L) endfor Reorder index set U to match the order of d(I_(u),L) for i = 1 to Ndo Select M images from I_(l),l ∈ L that are closest to I_(u) _(i)Register the selected atlases to I_(u) _(i) Generate a multi-atlassegmentation estimate of I_(u) _(i) end for Transfer the indices{u₁,...,u_(N)} from U to L end while

Multi-Atlas Propagation and Segmentation Refinement

Each label propagation is carried out by applying a modified version ofthe method for hippocampus segmentation described in van der Lijn et al.(2008). In this method the segmentations f^(j), j=1, . . . , M obtainedfrom registering M atlases are not fused to hard segmentation as inHeckemann et al. (2006) but are instead used to form a probabilisticatlas in the coordinate system of the target image I. This is an exampleof a “classifier fusion” technique.

In the original work, this subject-specific atlas is combined withpreviously learned intensity models for foreground and background togive an energy function that is optimized by graph cuts. We previouslyextended this method in a way that directly estimates the intensitymodels from the unseen image and that generalizes the approach to morethan one structure (Wolz et al., 2009). A Gaussian distribution for aparticular structure is estimated from all voxels which at least 95% ofthe atlases assign to this particular structure. The backgrounddistribution for a particular structure i with label f_(i) is estimatedfrom the Gaussian intensity distributions of all other structures withlabel f_(j), j≠i and of Gaussian distributions for the tissue classesT_(k), k=1, . . . , 3 in areas where no particular structure is defined.

By incorporating intensity information from the unseen image into thesegmentation process, errors obtained with conventional multi-atlassegmentation can be overcome.

Each registration used to build the subject-specific probabilistic atlasmay be carried out in three steps: rigid, affine and non-rigid. Rigidand affine registrations are carried out to correct for globaldifferences between the images. In the third step, two images arenon-rigidly aligned using a freeform deformation model in which aregular lattice of control point vectors are weighted using B-splinebasis functions to provide displacements at each location in the image(Rueckert et al., 1999). The deformation is driven by the normalizedmutual information (Studholme et al., 1999) of the pair of images. Thespacing of B-spline control points defines the local exibility of thenon-rigid registration. A sequence of control point spacings may be usedin a multi-resolution fashion (20 mm, 10 mm, 5 mm and 2.5 mm).

It will be appreciated that, in our method, we use multi-atlassegmentation to systematically label intermediate atlases that are thenused for multi-atlas segmentation of target images that are selectedaccording to their similarity with the previously labelled atlas images.Compared to previous work, we are dealing with a very diverse set ofimages. In such a scenario the gain from only registering similar imagesis more likely to outweigh the accumulation of registration errors.

Experimental Validation

We validated our new method experimentally as follows: We began bytaking an initial set of manually labelled atlases consisting of 30 MRimages from young and healthy subjects (age range 20-54, median age 30.5years) together with manual label maps defining 83 anatomical structuresof interest. In this set, the T1-weighted MR images had been acquiredwith a GE MR-scanner using an inversion recovery prepared fast spoiledgradient recall sequence with the following parameters: TE/TR 4.2 ms(fat and water in phase)/15.5 ms, time of inversion (TI) 450 ms, flipangle 20°, to obtain 124 slices of 1.5-mm thickness with a field of viewof 18×24 cm with a 192×256 image matrix.

We then used our method to propagate this initial set of atlases to adataset of 796 MR images acquired from patients with Alzheimer's Disease(AD) and mild cognitive impairment (MCI) as well as age matched controlsfrom the Alzheimer's Disease Neuroimaging Initiative (ADNI) database(www.loni.ucla.edu/ADNI). In the ADNI study, brain MR images had beenacquired at baseline and regular intervals from approximately 200cognitively normal older subjects, 400 subjects with MCI, and 200subjects with early AD.

From the results discussed below, it will be seen that this approachprovides more accurate segmentations due, at least in part, to theassociated reductions in inter-subject registration error.

Coordinate System Embedding

We applied the method for coordinate system embedding described above toa set of images containing the 30 initial atlases and the 796 ADNIimages. We used the first two features from spectral graph analysis toembed all images into a 2D coordinate system. The results of coordinatesystem embedding are displayed in FIG. 4. The original atlases form adistinct cluster on the left hand side of the graph at low values forthe first feature. Furthermore it can be seen that control subjects aremainly positioned at lower values, whereas the majority of AD subjectsis positioned at higher values. The hippocampal area for chosen examplesubjects is displayed in FIG. 4. These types of observations support theimpression that neighbourhoods in the coordinate system embeddingrepresent images that are similar in terms of hippocampal appearance.

All 796 images were segmented using five different approaches:

-   -   I. Direct segmentation using standard multi-atlas segmentation.    -   II. Direct segmentation using multi-atlas segmentation in        combination with an intensity refinement based on graph cuts.    -   III. Our new method, with M=30 and N=300 and no intensity        refinement after multiatlas segmentation.    -   IV. Our new method, with M=30 and N=1.    -   V. Our new method, with M=30 and N=300.

Evaluation of Segmentations

For evaluation we compared the automatic segmentation of the ADNI imageswith a manual hippocampus segmentation. This comparison was carried outfor all of the images for which ADNI provides a manual segmentation (182out of 796). Comparing these 182 subjects (Table 1) with the entirepopulation of 796 subjects (Table 2) shows that the subgroup ischaracteristic of the entire population in terms of age, sex, MMSE andpathology.

TABLE 1 Characteristics of the subjects used for comparison betweenmanual and automatic segmentation N M/F Age MMSE Normal 57 27/30 77.10 ±4.60 [70-89] 29.29 ± 0.76 [26-30] MCI 84 66/18 76.05 ± 6.77 [60-89]27.29 ± 3.22 [24-30] AD 41 21/20  76.08 ± 12.80 [57-88] 23.12 ± 1.79[20-26]

TABLE 2 Information relating to the subjects whose images were used inthis work N M/F Age MMSE Normal 222 106/216 76.00 ± 5.08 [60-90] 29.11 ±0.99 [25-30] MCI 392 138/254 74.68 ± 7.39 [55-90] 27.02 ± 1.79 [23-30]AD 182 91/91 75.84 ± 7.63 [55-91] 23.35 ± 2.00 [18-27]

An example for the segmentation of the right hippocampus of an ADsubject is shown in FIG. 5, with images (b), (c), (d), (e) correspondingto methods I, II, III and V respectively. A clear over-segmentation intoCSF space and especially an under-segmentation in the anterior part ofthe hippocampus can be observed, both in the case of multi-atlassegmentation with and without intensity-based refinement (methods I andII). The fact that the intensity-based refinement cannot compensate forthis error is due to the high spatial prior in this area that is causedby a significant misalignment of the majority of atlases in this area.The resulting high spatial prior cannot be overcome by theintensity-based correction scheme. When using the proposed frameworkwithout intensity-refinement (method III), the topological errors can beavoided, but the over-segmentation into CSF space is still present. Thefigure also shows that all observed problems can be avoided by using theproposed framework. In FIG. 5 (and also FIG. 6 and Table 3 below), theresults obtained using our new method are identified by the term “LEAP”(short for “Learning Embeddings for Atlas Propagation”).

The average overlaps as measured by the Dice coefficient or similarityindex (SI) (Dice, 1945) for the segmentation of left and righthippocampus on the 182 images used for evaluation are shown in Table 3.The difference between all pairs of the five methods is statisticallysignificant with p<0.001 on Student's two-tailed paired t-test.

TABLE 3 Dice overlaps for hippocampus segmentation Left hippocampusRight hippocampus Direct 0.775 ± 0.087 0.790 ± 0.080 [0.470-0.904][0.440-0.900] Direct, GC 0.820 ± 0.064 0.825 ± 0.065 [0.461-0.903][0.477-0.901] LEAP, N = 300, no GC 0.808 ± 0.054 0.814 ± 0.053[0.626-0.904] [0.626-0.900] LEAP, N = 1 0.838 ± 0.023 0.830 ± 0.024[0.774-0.888] [0.753-0.882] LEAP, N = 300 0.848 ± 0.033 0.848 ± 0.030[0.676-0.903] [0.729-0.905]

These results clearly show an improved segmentation accuracy androbustness for the proposed method. Our hypothesis is that by avoidingthe direct registration of images whose distance in the embedded spaceis too large but instead registering the images via multipleintermediate images improves significantly the segmentation accuracy androbustness of multi-atlas segmentation. To test this hypothesis we haveinvestigated the development of the segmentation accuracy as a functionof distances in the coordinate system embedding as well as the number ofintermediate steps. FIG. 6 shows this for the five segmentation methodsin the form of ten bar plots: Each bar plot corresponds to the averageSI overlap of 18 images (20 in the last plot). The first plot representsthe 18 images closest to the original atlases, the next plot representsimages slightly further from the original atlases and so on. Theseresults show the superiority of the proposed method over directmulti-atlas segmentation approaches in segmenting images that aredifferent from the original atlas set.

With increasing distance from the original atlases in the learnedmanifold, the accuracy of direct multi-atlas segmentation (method I) aswell as multiatlas segmentation with intensity-based refinement (methodII) steadily decreases. By contrast, our new method with both parametersettings shows a steady level of segmentation accuracy. It isinteresting to see that our method with a step width of N=1 (method IV)leads to worse results than the direct multiatlas methods up to acertain distance from the original atlases. This can be explained byregistration errors accumulated through many registration steps. Withincreasing distance from the atlases, however, the gain from usingintermediate templates, outweighs this registration error. Furthermore,the accumulated registration errors do not seem to increase dramaticallyafter a certain number of registrations. This is partly due to theintensity-based correction in every multi-atlas segmentation step whichcorrects for small registration errors. Segmenting the 300 closestimages with our new method before doing the next intermediate step(N=300, method V), leads to results at least as good as and often betterthan those given by the direct methods for images at all distances fromthe initial atlases. The importance of an intensity-based refinementstep after multi-atlas segmentation is also underlined by the results ofmethod III. When applying our new method without this step, the gaincompared to method I gets more and more significant with moreintermediate steps, but the accuracy still declines significantly whichcan be explained by a deterioration of the propagated atlases (note thatfor the first 300 images, method II and method V are identical, as aremethods I and III). The influence of N on the segmentation accuracy isgoverned by the trade-off between using atlases that are as close aspossible to the target image (small N) and using a design where aminimum number of intermediate steps are used to avoid the accumulationof registration errors (large N). Due to the computational complexity ofevaluating the framework, we restricted the evaluation to two values.

Volume Measurements

A reduction in hippocampal volume is a well-known factor associated withcognitive impairment (e.g. Jack et al. (1999); Reiman et al. (1998)). Tomeasure the ability of our method to discriminate clinical groups byhippocampal volume, we compared the volumes measured on the 182 manuallylabelled images to the ones obtained from our automatic method (methodV, LEAP with M=30 and N=300). Boxplots showing these volumes for theleft and right hippocampus are presented in FIG. 7, which shows averagehippocampal volumes for manual and automatic segmentation using methodIV. The discriminative power for the volume of left and righthippocampus between all pairs of clinical groups is statisticallysignificant with p<0.05 on a Student's t-test but is slightly lesssignificant than the manual discrimination.

FIG. 8 is a Bland-Altman plot showing the agreement between volumemeasurement based on manual and automatic segmentation of thehippocampus (method IV), with the solid line representing the mean andthe dashed lines representing ±1.96 standard deviations. This plotsupports the impression of the volume measures in FIG. 7 that theautomated method tends to slightly overestimate the hippocampal volumes.This over-segmentation is more significant for small hippocampi. Thesame phenomenon has been described for an automatic segmentation methodbefore by Hammers et al. (2007). The intraclass correlation coefficient(ICC) between the volume measurements based on the manual and automaticsegmentation is 0:898 (ICC (3,1) Shrout-Fleiss reliability (Shrout andFleiss, 1979)). This value is comparable to the value of 0:929 reportedin Niemann et al. (2000) for inter-rater reliability.

Discussion and Conclusion

In this work we have described our new method for propagating an initialset of brain atlases to a diverse population of unseen images viamultiatlas segmentation. We begin by embedding all atlas and targetimages in a coordinate system where similar images according to a chosenmeasure are close. The initial set of atlases is then propagated inseveral steps through the manifold represented by this coordinatesystem. This avoids the need to estimate large deformations betweenimages with significantly different anatomy and the correspondencebetween them is broken down into a sequence of comparatively smalldeformations. The formulation of the framework is general and is nottied to a particular similarity measure, coordinate embedding orregistration algorithm.

We applied our new method to a target dataset of 796 images acquiredfrom elderly dementia patients and age matched controls using a set of30 atlases of healthy young subjects. In this first application of themethod, we have applied it to the task of hippocampal segmentation eventhough the proposed framework can be applied to other anatomicalstructures as well. The proposed method shows consistently improvedsegmentation results compared to standard multi-atlas segmentation. Wehave also demonstrated a consistent level of accuracy for the proposedapproach with increasing distance from the initial set of atlases andtherefore with more intermediate registration steps. The accuracy ofstandard multi-atlas segmentation, on the other hand, steadilydecreases. This observation suggests three main conclusions:

(1) The decreasing accuracy of the standard multi-atlas segmentationsuggests that the coordinate system embedding used is meaningful. Theinitial atlases get less and less suitable for segmentation withincreasing distance.

(2) The almost constant accuracy of the proposed method suggests that,by using several small deformations, it is possible to indirectly deforman atlas appropriately to a target in a way that is not matched by adirect deformation within the multi-atlas segmentation framework used.

(3) The gain from restricting registrations to similar images countersthe accumulation of errors when using successive small deformations.

Our results indicate that, if many intermediate registrations are used,the segmentation accuracy initially declines quickly but then remainsrelatively constant with increasing distance from the initial atlases.The initial decline can be explained by an accumulation of registrationerrors which results from many intermediate registration steps. Thereason why the accuracy does not monotonically decline is likely to bedue to the incorporation of the intensity model during each multi-atlassegmentation step. By automatically correcting the propagatedsegmentation based on the image intensities, the quality of the atlascan be preserved to a certain level.

Apart from the obvious application of segmenting a dataset of diverseimages with a set of atlases based on a sub-population, the proposedmethod can be seen as an automatic method for generating a largerepository of atlases for subsequent multi-atlas segmentation with atlasselection (Aljabar et al., 2009). Since the manual generation of largeatlas databases is expensive, time-consuming and in many casesunfeasible, the proposed method could potentially be used toautomatically generate such a database.

Notwithstanding the challenge represented by variability due to imageacquisition protocols and inter-subject variability in a dataset aslarge and as diverse as the one in the ADNI-study, the results achievedwith our method compare well to state of the art methods applied to morerestricted datasets (van der Lijn et al., 2008; Morra et al., 2008;Chupin et al., 2009; Hammers et al., 2007) in terms of accuracy androbustness.

Summary

We have presented a new framework for the automatic propagation of a setof manually labelled brain atlases to a diverse set of images of apopulation of subjects. A manifold is learned from a coordinate systemembedding that allows the identification of neighbourhoods which containimages that are similar based on a chosen criterion. Within the newcoordinate system, the initial set of atlases is propagated to allimages through a succession of multi-atlas segmentation steps. Thisbreaks the problem of registering images which are very “dissimilar”down into a problem of registering a series of images which are“similar”. At the same time it allows the potentially large deformationbetween the images to be modelled as a sequence of several smallerdeformations.

Acknowledgement

The work leading to this invention has received funding from theEuropean Community's Seventh Framework Programme (FP7/2007-2011) undergrant agreement no. 224328.

REFERENCES

Aljabar, P., Heckemann, R., Hammers, A., Hajnal, J., Rueckert, D., 2009.Multi-atlas based segmentation of brain images: Atlas selection and itseffect on accuracy. NeuroImage 46 (3), 726-738.

Avants, B., Gee, J. C., 2004. Geodesic estimation for large deformationanatomical shape averaging and interpolation. NeuroImage 23 (Supplement1), S139-S150, Mathematics in Brain Imaging.

Bajcsy, R., Lieberson, R., Reivich, M., August 1983. A computerizedsystem for the elastic matching of deformed radiographic images toidealized atlas images. J. Comput. Assisted Tomogr. 7, 618-625.

Blezek, D. J., Miller, J. V., 2007. Atlas stratification. Medical ImageAnalysis 11 (5), 443-457.

Chung, F. R. K., 1997. Spectral graph theory. Regional Conference Seriesin Mathematics, American Mathematical Society 92, 1-212.

Chupin, M., Hammers, A., Liu, R., Colliot, O., Burdett, J., Bardinet,E., Duncan, J., Garnero, L., Lemieux, L., 2009. Automatic segmentationof the hippocampus and the amygdala driven by hybrid constraints: Methodand validation. NeuroImage 46 (3), 749-761.

Collins, D. L., Holmes, C. J., Peters, T. M., Evans, A. C., 1995.Automatic 3-D model-based neuroanatomical segmentation. Human BrainMapping 3 (3), 190-208.

Davis, B. C., Fletcher, P. T., Bullitt, E., Joshi, S., 2007. Populationshape regression from random design data. In: ICCV. pp. 1-7.

Dice, L. R., 1945. Measures of the amount of ecologic associationbetween species. Ecology 26 (3), 297-302.

Ericsson, A., Aljabar, P., Rueckert, D., 2008. Construction of apatient-specific atlas of the brain: Application to normal aging. In:ISBI. IEEE, pp. 480-483.

Gee, J. C., Reivich, M., Bajcsy, R., March-April 1993. Elasticallydeforming 3D atlas to match anatomical brain images. Journal of ComputerAssisted Tomography 17 (2), 225-236.

Hammers, A., Allom, R., Koepp et al., M. J., August 2003.Three-dimensional maximum probability atlas of the human brain, withparticular reference to the temporal lobe. Human Brain Mapping 19 (4),224-247.

Hammers, A., Heckemann, R., Koepp, M. J., Duncan, J. S., Hajnal, J. V.,Rueckert, D., Aljabar, P., 2007. Automatic detection and quantificationof hippocampal atrophy on MRI in temporal lobe epilepsy: Aproof-of-principle study. NeuroImage 36 (1), 38-47.

Heckemann, R. A., Hajnal, J. V., Aljabar, P., Rueckert, D., Hammers, A.,2006. Automatic anatomical brain MRI segmentation combining labelpropagation and decision fusion. NeuroImage 33 (1), 115-126.

Jack, C. R., J., Petersen, R. C., Xu, Y. C., O'Brien, P. C., Smith, G.E., Ivnik, R. J., Boeve, B. F., Waring, S. C., Tangalos, E. G., Kokmen,E., 1999. Prediction of AD with MRI-based hippocampal volume in mildcognitive impairment. Neurology 52 (7), 1397-1407.

Joshi, S., Davis, B., Jomier, M., Gerig, G., 2004. Unbiaseddiffeomorphic atlas construction for computational anatomy. NeuroImage23 (Supplement 1), 151-160, mathematics in Brain Imaging.

Mazziotta, J. C., Toga, A. W., Evans, A. C., Fox, P. T., Lancaster, J.L., June 1995. A probabilistic atlas of the human brain: Theory andrationale for its development. the international consortium for brainmapping (ICBM). NeuroImage 2 (2a), 89-101.

Morra, J. H., Tu, Z., Apostolova, L. G., Green, A. E., Avedissian, C.,Madsen, S. K., Parikshak, N., Hua, X., Toga, A. W., Jr., C. R. J.,Weiner, M. W., Thompson, P. M., 2008. Validation of a fully automated 3Dhippocampal segmentation method using subjects with Alzheimer's diseasemild cognitive impairment, and elderly controls. NeuroImage 43 (1),59-68.

Niemann, K., Hammers, A., Coenen, V. A., Thron, A., Klosterktter, J.,2000. Evidence of a smaller left hippocampus and left temporal horn inboth patients with first episode schizophrenia and normal controlsubjects. Psychiatry Research: Neuroimaging 99 (2), 93-110.

Reiman, E. M., Uecker, A., Caselli, R. J., Lewis, S., Bandy, D., deLeon, M. J., Santi, S. D., Convit, A., Osborne, D., Weaver, A.,Thibodeau, S. N., August 1998. Hippocampal volumes in cognitively normalpersons at genetic risk for Alzheimer's disease. Annals of Neurology 44(2), 288-291.

Rohlfing, T., Russakoff, D. B., Jr, C. R. M., 2004. Performance-basedclassifier combination in atlas-based image segmentation usingexpectation-maximization parameter estimation. IEEE Trans. Med. Imaging23 (8), 983-994.

Rueckert, D., Sonoda, L. I., Hayes, C., Hill, D. L. G., Leach, M. O.,Hawkes, D. J., August 1999. Nonrigid registration using free-formdeformations: Application to breast MR images. IEEE Trans. MedicalImaging 18 (8), 712-721.

Sabuncu, M. R., Balci, S. K., Golland, P., 2008. Discovering modes of animage population through mixture modeling. In: MICCAI (2). Vol. 5242 ofLecture Notes in Computer Science. Springer, pp. 381-389.

Shrout, P., Fleiss, J., 1979. Intraclass correlation: uses in assessingrater reliability. Psychological Bulletin 86, 420-428.

Studholme, C., Hill, D. L. G., Hawkes, D. J., January 1999. An overlapinvariant entropy measure of 3D medical image alignment. PatternRecognition 32 (1), 71-86.

Tang, S., Fan, Y., Kim, M., Shen, D., February 2009. RABBIT: rapidalignment of brains by building intermediate templates. In: SPIE. Vol.7259.

van der Lijn, F., den Heijer, T., Breteler, M. M., Niessen, W. J., 2008.Hippocampus segmentation in MR images using atlas registration, voxelclassification, and graph cuts. NeuroImage 43 (4), 708-720.

van Rikxoort, E. M., Isgum, I., Staring, M., Klein, S., van Ginneken,B., 2008. Adaptive local multi-atlas segmentation: application to heartsegmentation in chest CT scans. Vol. 6914. SPIE, p. 691407.

von Luxburg, U., 2007. A tutorial on spectral clustering. Statistics andComputing 17 (4), 395-416.

Warfield, S. K., Zou, K. H., Wells III, W. M., 2004. Simultaneous truthand performance level estimation (STAPLE): an algorithm for thevalidation of image segmentation. IEEE Trans. Med. Imaging 23 (7),903-921.

Wolz, R., Aljabar, P., Heckemann, R. A., Hammers, A., Rueckert, D.,2009. Segmentation of subcortical structures and the hippocampus inbrain MRI using graph-cuts and subject-specific a-priori information.IEEE International Symposium on Biomedical Imaging—ISBI 2009.

Wu, M., Rosano, C., Lopez-Garcia, P., Carter, C. S., Aizenstein, H. J.,2007. Optimum template selection for atlas-based segmentation.NeuroImage 34 (4), 1612-1618.

1-26. (canceled)
 27. A method of processing medical images, performed bya computer processor and comprising the steps of: (a) obtaining one ormore atlases containing one or more images in which one or moreanatomical features have been labelled with label data; (b) obtaining aplurality of unlabelled images; (c) comparing the labelled andunlabelled images and selecting one or more unlabelled images that mostclosely resemble(s) one or more of the labelled images; (d) to each ofthose selected image(s), propagating label data from one or more of theclosest of the labelled images, thereby labelling the correspondinganatomical feature(s) of each of the selected image(s) and causing theselected image(s) to become labelled image(s); and (e) iterativelyrepeating from step (c), thereby labelling others of the unlabelledimages.
 28. A method as claimed in claim 27, wherein the step ofcomparing the labelled and unlabelled images comprises embedding theimages into a low-dimensional coordinate system.
 29. A method as claimedin claim 27, wherein the low-dimensional coordinate system is atwo-dimensional coordinate space.
 30. A method as claimed in claim 27,wherein the step of comparing the labelled and unlabelled imagescomprises defining a set of pairwise measures of similarity by comparingone or more respective anatomical features for each pair of images inthe set of images.
 31. A method as claimed in claim 30, wherein the stepof comparing the labelled and unlabelled images further comprisesperforming a spectral analysis operation on the pairwise measures ofsimilarity.
 32. A method as claimed in claim 30, wherein the pairwisemeasures of similarity represent the intensity similarity between a pairof images.
 33. A method as claimed in claim 30, wherein the pairwisemeasures of similarity represent the amount of deformation between apair of images.
 34. A method as claimed in claim 27, wherein the step ofpropagating label data comprises propagating label data from a pluralityof the closest of the labelled images, based on a classifier fusiontechnique.
 35. A method as claimed in claim 27, further comprising,after step (d) and before step (e), a step of performing anintensity-based refinement operation on the newly-propagated label data.36. A method as claimed in claim 27, wherein the images are of differentsubjects.
 37. A method as claimed in claim 27, wherein at least some ofthe images are of the same subject but taken at different points intime.
 38. A method as claimed in claim 27, wherein the images aremagnetic resonance images.
 39. A method as claimed in claim 27, furthercomprising labelling an anatomical feature representative of thepresence or absence of a condition and using that feature to derive abiomarker for that condition.
 40. A method as claimed in claim 39,further comprising allocating a subject to a diagnostic category on thebasis of the biomarker.
 41. A method as claimed in claim 39, furthercomprising quantifying a subject's response to treatment on the basis ofthe biomarker.
 42. A method as claimed in claim 39, further comprisingselecting a subject's treatment on the basis of the biomarker.
 43. Acomputer system, image processing apparatus, or imaging apparatusarranged to implement a method as claimed in claim
 27. 44. An imagingapparatus as claimed in claim 43, being a medical scanner.
 45. Animaging apparatus as claimed in claim 44, being an MRI scanner.
 46. Acomputer program comprising coded instructions for implementing a methodas claimed in claim 27 when run on a computer processor, optionallywherein the computer program is encoded in a computer-readable medium orphysical carrier signal.